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Nobuyuki Kanematsu 

Department of Accelerator and Medical Physics, National Institute of Radiological Sciences, 4-9-1 Anagawa, Inage-ku, Chiba 263-8555, Japan 



oo 

o 

O 

(N 
•*-> 
O 

O 

oo 



Abstract 

This study provides an accurate, efficient, and simple multiple scattering formulation for heavy charged particles such as protons 
and heavier ions with a new form of scattering power that is a key quantity for beam transport in matter. The Highland formula 
for multiple scattering angle was modified to a scattering-power formula to be used within the Fermi-Eyges theory in the presence 
of heterogeneity. An analytical formula for RMS end-point displacement in homogeneous matter was also derived for arbitrary 
ions. The formulation was examined in terms of RMS angles and displacements in comparison with other formulations and 
measurements. The results for protons, helium ions, and carbon ions in water agreed with them at a level of 2% or the differences 
were discussed. 
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Theory of Coulomb scattering was developed in detail many 
years ago. Rutherford [1] studied elastic scattering of a charged 
particle by a point-like nucleus and found 1/0 4 behavior in the 
distribution of scattering angle 9. In matter, a particle under- 
goes multiple processes of scattering by nuclei with electric 
field screened by orbital electrons. Moliere developed a rig- 
orous theory for such a system and formulated an analytical 
expression of the angular distribution for a particle interacting 
with target atoms [2]. The resultant Moliere distribution has 
Gaussian behavior at small angles and 1/0 4 behavior at large 
angles. The Gaussian behavior is formed by multiple small- 
angle scatters with the central-limit theorem in statistics, while 
the I/O 4 tail reflects single large-angle scattering. 

In cases of accelerated particles, either whole or subset of 
them are often modeled as a Gaussian beam because it is phys- 
ically approximate to the reality, algorithmically efficient for 
good localization nature, and numerically easy with the stan- 
dard math library. Fermi and Eyges [3] developed a theory for 
Gaussian beam transport in matter, where number density <t> at 
transverse position y and angle 9 for a -particle system is de- 
scribed as 
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which is characterized by angular variance 9 2 — Yji=\ 9] IN, spa- 
tial variance y 2 = J^f-iyf/N, and covariance y0 - Y!i=\y$ilN ■ 
As the particles undergo energy loss and multiple scattering, 
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where scattering power T is the key quantity that drives the 
beam development in heterogeneous system along longitudinal 
position x. Note that we consistently deal with projected posi- 
tion y and angle 6 in this work and that 9 2 , y9, y 2 , T, and con- 
stant E 2 (Sec. 2.1) are thus 1/2 of the conventional definitions 
with radial positions and polar angles. We also give specific 
symbols to such quantities later in various formulations. 

For a given model of 9 2 (x), it is possible to numerically ob- 
tain effective scattering power f = [9 2 (x + Ax) - 9 2 (x)]/Ax for 
small step Ax [4, 5]. However, it is desirable to have an analyt- 
ical formula for the scattering power, which is the main objec- 
tive of this work, not only for theoretical cleanness but also for 
further analytical derivation of physical quantities. 

In the following sections, we examine various Gaussian ap- 
proximation models for multiple scattering, propose a new form 
of scattering power as a better solution, and compare with other 
studies in terms of resultant physical quantities such as angles 
and displacements. Although this work is primarily intended 
for heavy charged particle radiotherapy with protons or heavier 
ions, the subject is substantially general and may be useful for 
other beam-transport applications [6, 7]. 
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2. Materials and methods 

2.1. Gaussian approximations for multiple scattering 

Fermi-Rossi (FR) formulation. Fermi and Rossi [8] developed 
a theory of multiple scattering with RMS angle #fr and scatter- 
ing power Tfr formulated as 
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where E s = m e c 2 y/2n:/a ~ 15.0 MeV is a constant energy, z, 
p, and v are the charge/e, the momentum, and the speed of the 
particle, and Xq is the radiation length that conveniently encap- 
sulates the material properties [9]. Since the F-R formulation 
relies on the central-limit theorem, existence of the I/O 4 tail 
at large angles leads to inaccuracy. Nevertheless, the FR or 
equivalent form of T oc (V/ pv) 2 has been commonly used in the 
Fermi-Eyges theory [3, 10, 11]. 

Moliere-Hanson (MH) formulation. In the Moliere theory, 
the angle distribution has a Gaussian term with RMS angle 
(Xc/ V2) Vb, where Xc an d B can be interpreted as the charac- 
teristic angle per scattering and the mean number of scattering 
per particle. The width of the central part of the Moliere distri- 
bution is, however, slightly narrower than that of the Gaussian 
term due to contributions of non-Gaussian terms. Hanson et al. 
[12] found the best-approximate Gaussian RMS angle for the 
central part as 

0MH = ^ Vfl-1.2. (6) 
V2 

Although the MH angle (9mh is often used in Gaussian beam 
models [4, 13, 14], the complexity of the theory would discour- 
age its direct use in demanding applications and it is difficult to 
formulate an analytical scattering power to handle heterogene- 
ity. As to the details of the Moliere theory including Fano cor- 
rection, this work strictly follows the formalism by Gottschalk 
etal. [15]. 

Integral Highland ( iH) formulation. Highland [ 16] introduced 
a simple correction term to the integral form of the FR formula 
with an optimized energy constant, 1 for better agreement with 
the MH angle. Gottschalk et al. [15] then generalized for thick 
targets with formula 
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where lg = log 10 is the common logarithmic function and x is 
the thickness of a homogeneous target. The iH angle 6>;h was 
experimentally verified to be accurate [15, 17] and has been 
used in practice [18, 19]. 

Often, a target may have composite structure of multiple 
elements. In such cases, use of the quadratic additivity rule 



(9 2 H = 2i ma Y appear natural, but is incompatible with the 
ill-behaved logarithmic term for thin target layers to handle het- 
erogeneity. Kanematsu et al. [5] addressed the problem by fur- 
ther generalization 
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where x is interpreted as the longitudinal position in the target 
and 9m is the RMS angle growing with x. Radiative path length 
t is defined as 

^to= I (9) 
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for any composite target or heterogeneous system with radia- 
tion length Xo varying with x. However, it is also difficult to de- 
rive an exact scattering power from Eq. (8) due to the involved 
integral terms. 

0verds-Schneider ( 0S) formulation. 0veras [20] found a good 
approximate relation with pv, 
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where R is the residual range of the particle expected in water, 
k — 1.08 (Sec. 2.3.1) is a constant, and po, vq, and Rq are the 
initial values on the incidence. While R and k were originally 
formulated as material-dependent, 2 we determined to take wa- 
ter as a reference material to measure kinetic energy E or pv. 

Schneider et al. [21] applied 0veras's relation to the FR for- 
mulation and proposed a RMS angle with correction to best 
reproduce Gottschalk's [15] experimental data, 
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-4.86, (12) 



for the target material with density p, radiation length Xo, and 
stopping -power ratio ps with respect to water. In the 0S formu- 
lation, the scattering power is analytically given by 
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The 0S formulation does not seem to be popular despite the 
potential adaptability to the Fermi-Eyges theory. 



Highland [16] quantified the constant in a different form as 17.5 MeV that 
would be 13.9 MeV in the standard form, whereas 14.1 MeV is commonly 
referred to as Highland's constant and has been used as the standard. 



Schneider et al. [21] found strong correlation between k and the radiation 
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Figure 1: The differential Highland (solid) and the 0veras-Schneider (dotted) 
scattering powers relative to the Fermi-Rossi scattering power (dashed) of water 
as a function of normalized depth for particles with incident range Rq = Xq . 



2.2. Alternative scattering power 

The essence of the Highland formulation is the presence of 
the correction factor to the integral form of the FR formula (5), 
which is considered as a variable-separation approximation. It 
is thus natural to introduce a similar correction factor /dH in the 
differential Highland (dH) formula or the dH scattering power 
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to accommodate the single-scattering effect. The average of 
the factor for the entire path should coincide with the original 
Highland correction squared or the ratio f?j H /#p R , 

1 f' , , / lgn 2 /l4.1MeV\ 2 

jjf /WO* + (-5-). (15) 

leading to a simple solution, 
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with which T^h is exactly defined. 

Figure 1 shows relative strength of scattering powers 7Vr, 
Tan, and 7gjs of water (p = p s = 1, Xq = 36.08 cm) as a 
function of normalized depth p^x/Ro = 1 - R/Rq for particles 
with incident range Ro = Xq. These curves exhibit significant 
differences reflecting the complexity of the formulations. 

2.3. Beam development 

2.3.1. Range-energy relation 

For a semi-relativistic heavy charged particle, the Bethe the- 
ory describes its mean stopping behavior in matter [22]. At 
fixed velocity v, particle dependency of stopping power S = 
-dE/dx oc zr and that of kinetic energy E oc mass m lead 



that of residual range to R — dE' /S(E') oc m/z 2 . Relation 
pv = E{E + 2mc 2 )/(E + mc 2 ) leads to pv oc m. These extend 
0veras's relation for fully stripped ions or nuclei with charge 
z e and mass m = A h as 
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where m p = 1.0073 u is the proton mass and u — 931.5 MeV/c 2 
is the atomic mass unit. Parameters k and A were deter- 
mined with ICRU [22] data points for protons in water with 
(£/MeV,fl/cm) = (200, 25.96) and (400,82.25). 

2.3.2. Numerical computation for heterogeneous systems 

The beam development is computed in a stepwise manner to 
deal with particle energy loss and medium heterogeneity. In 
small step Ax from x, the residual range and the radiative path 
length are modified by 



AR = -ps(x) Ax, 
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and integrals (2)-(4) are translated into increments 
Afi = f Ax, 
Ay6 ■■ 
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Effective scattering power T in the FR formulation is 
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where the geometric-mean pv represent its effective value ac- 
curately for steps within 20% of the residual range [15]. The 
effective 0S and dH scattering powers are 
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where R as R + AR/2 and I « I + A£/2 are approximate effective 
values for the step. 

Since the R-E relation may have 1 % or more uncertainty, we 
limit depth step ps Ax to be more than 0.5% of the initial range 
Rq to balance accuracy and efficiency. In addition, we limit 
the step to be less than 10% of the residual range R and within 
distance 5 to the next heterogeneous layer boundary, by 
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Table 1 : Atomic properties (mass density, mass-electron density, mean excita- 
tion energy, radiation mass length, and effective mass-stopping-power ratio) of 
water and target materials. 



Material 


Pl-h 

r ' cm- 1 


n e i mol 


7/eV 




Ps j cm 1 
P g 


Water 


l 


0.5551 


75 


36.08 


1 


Beryllium 


1.85 


0.4438 


63.7 


65.19 


0.8195 


Copper 


8.96 


0.4564 


322 


12.86 


0.6674 


Lead 


11.35 


0.3958 


823 


6.37 


0.4913 



In the last step Ax = R/ps to the end point, T would diverge 
as pv — > in Eq. (22). Since angle 9 loses its physical signif- 
icance there, we should only deal with y 2 . The last increment 
Ay 2 can be analytically calculated with Eqs. (4), (14), and (17), 
as 

(26) 
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2.3.3. Analytical formula for homogeneous systems 

Analytical integral (26) for the last step is valid for larger 
distances in the absence of heterogeneity. Suppose an infinites- 
imal parallel beam with initial range Rq is incident into a homo- 
geneous target, the particles traverse distance Ro/ps and stop 
with mean square displacement y 2 . The effective Highland 
correction factor to y 2 should be 7? 2 -weighted mean of fm({) 
in Eq. (16) with t = (R - R)/(p s X ) as 
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We thus obtain analytical RMS end-point displacement 
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for stopping ions in a homogeneous system. 



(28) 



2.4. Application and validation 

Range-energy relation. Since this work is heavily dependent 
on the R-E (or pv) relation, we first examined its accuracy 
against the standard data [22] for the interested energy re- 
gion of E/A < 400 MeV in comparison with R-E relation 
R/cm = 0.0022 (£/MeV) 177 similarly proposed by Bortfeld 
[23] for E < 250 MeV protons. 

Scattering angle. RMS angles of the various formulations 
were compared against the reference MH angles calculated by 
Gottschalk et al. [15] for 158.6 MeV protons incident into 
beryllium, copper, and lead targets with properties in Table 1, 
where the effective mass-stopping -power ratio ps/p of the ma- 
terials were derived from the 158.6 MeV proton ranges in the 
targets [15] and that in water Ro = 17.30 cm estimated by rela- 
tion (17). 



Transverse displacement. RMS transverse displacements <r y in 
water (ps = \,Xq- 36.08 cm) were calculated with the FR, 0S, 
and dH scattering powers as a function of depth x for projectile 
nuclei with incident range Rq = 29.4 cm 'H, 29.4 cm 4 He, and 
29.7 cm 12 C to compare with Phillips's measurements of 1/e 
radius V2cr v [11]. 

End-point displacement. Similarly, the RMS transverse dis- 
placements at the end point, cr y0 , of 'H, 4 He, and 12 C nuclei 
incident into water were calculated with varied incident ranges 
and were compared with the analytical formula (28) and mea- 
surements. 

In addition to Phillips's measurements [11], we included two 
proton points measured by Preston and Kohler in their unpub- 
lished work in 1968, which were V2cr y0 = (0.346 + 0.009) cm 
forflo = 11.4 cm and V2cr v = (0.368 + 0.010) cm at x = 12.4 
cm for Rq = 12.8 cm converted to V2cr v o = 0.391 cm with 
their universal curve 
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where xr = ps x/Ro is the normalized depth. We also added 
other two points calculated with the Hanson form of Moliere's 
theory by Deasy [13] for 160 and 250 MeV protons. Deasy 
evaluated transversal FWHM's of Bragg peak as 2 V21n2 <x v = 
0.91 cm at x — 17.5 cm for Rq = 17.65 cm and 2 V2 In 2 cr v = 
1.96 cm at unspecified depth for Ro = 37.94 cm. The conver- 
sion factor for the former by the universal curve is 1.017 that 
is also applied to the latter for the best guess, leading to <x v o = 
0.393 cm and 0.846 cm. 

Heterogeneity handling. We calculated behaviors of Ro = 29.4 
cm protons in a multilayered bi-density water target, where 
layer ; = {1,2,3,4,...} of common thickness t and density 
Pi = [1.1,0.9, 1.1,0.9, ...} was placed with its upstream face at 
Xj = {0, t, 2t, 3t, ...}. The RMS end-point displacements were 
calculated by numerical integral (21) using the FR, 0S, and dH 
scattering powers (5), (13), and (14) and the effective scattering 
power fin per layer derived from the Gottschalk form of RMS 
angle iH (7), 

t \ 9 B X 0w j Jo \pv( Xi + dt')) tX 0w 
The thickness t was varied in the extent of 0.01-1 cm. 



3. Results 

Range-energy relation. Figure 2 shows that the R-pv relation 
(17) and the standard ICRU data [22] for protons agreed within 
either 0. 1 cm or 1 % for 0-400 MeV, which would not have been 
accomplished with Bortfeld's relation. 

Scattering angle. The RMS scattering angles are compared in 
Fig. 3 and in Table 2. For the homogeneous target systems, 
the two variations of Highland angles f?;H and 9m were in fact 
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Figure 2: Range-energy relation curve for protons in water translated from the 
R—pv relation (solid), Bortfeld's curve (dashed), and ICRU data points (o). 
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Figure 3: Calculated RMS angles relative to the the Moliere-Hanson an- 
gles 0mh [15] (dashed) by Fermi-Rossi (x), integral Highland (+), 0veras- 
Schneider (A), and differential Highland (o) formulations for 158.6 MeV pro- 
tons scattered by (a) beryllium, (b) copper, and (c) lead targets as a function of 
normalized target thickness psx/Rg. 



Table 2: Calculated RMS angles for 158.6 MeV protons by beryllium, copper, 
and lead targets of normalized thicknesses psx/Rq = 1% and 10%, by Fermi- 
Rossi (#frX Moliere-Hanson (6mh) [15], integral Highland ($ih), 0veras- 
Schneider (00s)> and differential Highland (Ban) formulations. 



Target 
Psx/Ro 



Beryllium 



Copper 



Lead 



1% 10% 1% 10% 1% 10% 



(Wmrad 2.93 9.49 7.23 23.5 11.9 38.6 

fWmrad 2.01 7.17 5.63 20.4 9.75 35.8 

(Wmrad 1.99 7.45 5.51 20.3 9.60 35.2 

% s /mrad 2.02 7.06 5.53 20.0 9.26 33.6 

dH /mrad 2.04 7.61 5.63 20.7 9.78 35.8 




Figure 4: RMS transverse displacements in water for projectiles (a) Rq = 29 A 
cm 1 H, (b) 29.4 cm 4 He, and (c) 29.7 cm 12 C as a function of depth, and (d) all 
of them in normalized scale, by numerical computation with the Fermi-Rossi 
(dashed), differential Highland (solid), and 0veas-Schneider (dotted) scattering 
powers. Markers indicate Phillips's measurements. 



equally consistent with the MH angle while the FR angle de- 
viated from them by about 50% at small thicknesses. The 0S 
angle was superior for the beryllium target and inferior for the 
copper and lead targets to the Highland angles in terms of agree- 
ment with the MH angle. 

Transverse displacement. Figures 4 (a)-(c) show the growths 
of transverse displacement of : H, 4 He, and 12 C nuclei in water. 
In terms of relative agreement with the measurements, the dH 
and 0S formulations were excellent for 1 H, the FR and dH for- 
mulations were good for 4 He, and the FR was the excellent for 
12 C. Considering their absolute scale and inherent experimental 
difficulties, the best would be the dH formulation with agree- 
ment within 2% or 0.02 cm everywhere. Figure 4 (d) shows 
the behaviors in the self-normalized scale, indicating that the 
relative displacement is very insensitive to the scattering for- 
mulations. The universal curve (29) would coincide with these 
curves. 



End-point displacement. Figures 5 (a)-(c) show the RMS dis- 
placements of 'H, 4 He, and 12 C nuclei at the end point in water 
for varied incident ranges. In terms of agreement with the mea- 
surements, the dH formulation would be the best for the same 
reason as for o~ y drawn mostly with the same data. The ana- 
lytical dH curve agreed well with the numerical computation. 
As shown in Fig. 5 (d), behavior of the dH and 0S formula- 
tions were approximately linear with relation <x v0 0.023 Ro, 
while the FR formulation deviated by 10% or more with larger 
non-linearity. 
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Figure 5: RMS end-point displacements in water for (a) 'H, (b) 4 He, and (c) 
12 C, and (d) displacement-range ratio for 'H as a function of incident range Rq 
by numerical computation with the Fermi-Rossi (dashed), differential Highland 
(solid), and 0veas-Schneider (dotted) scattering powers, the analytical formula 
(thick gray), measurements by Phillips (X, o, +) and by Preston and Kohler (A), 
and Moliere-Hanson calculations by Deasy (o). 
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Figure 6: RMS end-point displacements of Rq = 29.4 cm protons in a tar- 
get comprised of alternative p = 1.1 and 0.9 water layers of varied thickness 
calculated with the Fermi-Rossi (x), 0veas-Schneider (A), and differential 
Highland (o) scattering powers, and the integral Highland angle per layer (+). 



Heterogeneity handling. Figure 6 shows the behaviors of the 
various formulations against heterogeneity size. The layer-wise 
calculation of the iH angle (30) actually caused large error as 
addressed in Sec. 2.1. At typical layer thickness t = 0.1 cm, 
the misused iH formulation resulted in 26% underestimation in 
RMS displacement with respect to the dH formulation. 

4. Discussion 

The Fermi-Eyges theory generally, smartly, and efficiently 
describes a Gaussian beam with a set of a few beam-defining 
parameters computed by path integrals. For radiotherapy, it is 
useful to deal with field formation in beam delivery systems, 
beam customization for individual treatment targets, and beam 
transport in patients to give a variable pencil kernel for dose 
convolution algorithms. [5, 24, 25]. The scattering power de- 
termines the accuracy of the modeled Gaussian beam. Fortu- 
nately, the relative error will be about a half in RMS angle and 
displacement because of relations T cc Q 1 cc y 1 . This may be 
one of the reasons why the apparent differences in T were not 
very significant in cr v in this work. 

The FR formulation (5) without consideration of the single- 
scattering effect was inaccurate for thin targets of beryllium, 
copper, and lead. The inaccuracy or the single-scattering effect 
decreased with the target thickness. This is analogous to the 
relative decrease of Hanson's correction with increasing B in 
Eq. (6) and could be explained as follows. Increase of the multi- 
ple scattering with thickness also increases the threshold to dis- 
tinguish large-angle scattering, and thus decreases the single- 
scattering effect. For particles stopping in water, the inaccuracy 
was in fact small. 

The Highland formula is an empirical approximation of Han- 
son's another empirical approximation of the Moliere theory. 
There are even variations of constant parameters within the for- 
malism [26]. Nevertheless, the standard form adopted in this 
work has been experimentally validated [15, 17, 18, 27, 28] 
and thus should be a reasonable choice. The original Highland 
formula, however, was designed for homogeneous systems and 
causes large errors when misused for fine heterogeneity. 

We formulated a scattering power (14) by approximate dif- 
ferentiation of the generalized Highland formula (8). The two 
forms were verified to be equivalent. We analytically derived 
a general formula (28) for RMS end-point displacement cr v0 of 
stopping ions in a homogeneous system, which showed surpris- 
ingly linear behavior with incident range Rq. Since the formula 
has separate factors of Rq, z, and m dependencies, one can easily 
compare the scattering effects for different ion beams. 

The R-pv relation in the 0S formulation is superior to the 
conventional R-E relation and was fully utilized in this work. 
As to the scattering, the 0S formulation generally resulted in 
smaller displacements in water than those of the Highland for- 
mulation by several percent, as so found by Schneider et al. 
[21]. In terms of agreement with measurements by Phillips, 
measurements by Preston and Kohler, and MH calculations by 
Deasy, the 0S formulation was slightly inferior to the Highland 
formulation, although it was a fit to the other experiment. 
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While both formulations deal with the effect of large-angle 
scattering in Gaussian approximation, there is an essential dif- 
ference. Highland's correction is based on radiative path length 
I — L dx' /Xq(x') that may represent the multiple scattering ac- 
cumulated for the entire path. This approach sounds reasonable 
because large angles can only be defined with respect to the 
multiple-scattering angle. In contrast, Schneider's correction 
is based on relative residual range R/Rq and material-specific 
constants cq and c\, which are instantaneous quantities. That 
approach seems to assume invariance of atomic composition 
and could be inappropriate for systems with multiple materials 
of different compositions. 

5. Conclusions 

We formulated a scattering power with correction for single- 
scattering effect based on the Highland formula. It can be gen- 
erally used within the Fermi-Eyges theory for beam applica- 
tions with heavy charged particles in the presence of hetero- 
geneity. We derived an analytical formula for RMS displace- 
ment of ions stopped in a homogeneous system, which showed 
very linear behavior with the incident range in water. 

The single-scattering effect was as large as 50% for thin tar- 
gets and generally decreases with thickness. For ions stopping 
in water, the effect in displacement was about 10% and the rela- 
tive growth with depth is generally insensitive to the scattering- 
power formulations. The numerical and analytical calculations 
of this work were consistent and agreed with other experimental 
studies within 2% or 0.02 cm in RMS displacement. 
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